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This paper describes NASA Lewis' studies of complex reacting systems at 
high temperature. The changes which occur are the result of many different 
chemical reactions occurring at the same time. Both an experimental and a 
theoretical approach are needed to fully understand what happens In these sys- 
tems. The latter approach Is discussed herein. We present the differential 
equations which describe the chemical and thermodynamic changes, and we 
describe their solution by numerical techniques using a detailed chemical 
mechanism. Several different comparisons of computed results with experimen- 
tal measurements are also given. These Include the computation of (1) species 
concentration profiles In batch and flow reactions, (2) rocket performance In 
nozzle expansions, and (3) pressure versus time profiles In hydrocarbon Igni- 
tion processes. The examples Illustrate the use of detailed kinetic computa- 
tions to elucidate a chemical mechanism and to compute practical quantities 
such as rocket performance. Ignition delay times, and Ignition lengths In flow 
processes. 


INTRODUCTION 

For many years Lewis has been studying complex, reacting gas-phase sys- 
tems at high temperature. The changes that occur In such systems are the 
result of many Individual chemical reactions occurring at the same time. We 
use a two-step approach to understand what happens In such a system: (1) 

experimental measurement of the temporal changes In temperature, pressure, or 
composition of the system, and (2) theoretical computation of these changes In 
an attempt to match the experimental results. This report describes our work 
In the computation of complex-system chemical kinetics as focused on combus- 
tion systems. Brabbs describes our experimental efforts In his paper for this 
symposium. Early attempts to compute the progress of complex reactions 
Involved various simplifications. These were of two types: (1) one or two 
global reactions were substituted for the actual set of many Individual, 
simultaneous reactions which occur In a gaseous system with several species 
present, and (2) reactions Involving some of the very reactive atoms and radi- 
cals were assumed to be very fast and, therefore, always In chemical equilib- 
rium. Differential equations for the rates of change of a few key species 
were solved analytically. Although both of these approaches have had limited 
success, they only explain some of the general features of a complex reaction. 
They are, at best, simplified approximations of the actual process, which are 
only valid under a limited set of conditions. The most fruitful approach to 
complex kinetics computations Is to numerically Integrate a system of differ- 
ential equations derived from the laws of conservation of mass, energy, and 
(If flow Is Involved) momentum. These differential equations Involve the 
rates of change of species concentrations oj (moles per unit mass of mix- 
ture) and temperature T, and they may also Involve density p and velocity 
V for a flow process. When these differential equations are solved, a set of 
Individual reactions Is assumed to occur simultaneously among all the species, 
and these reactions are usually assumed to be reversible (with the ratio of 
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the forward to reverse rate constant being equal to the equilibrium constant 
for any reaction) . 

There are two ways to use these detailed kinetic computations. First, 
they may be used to elucidate a reaction mechanism or to determine a single 
unknown rate constant In conjunction with experimental data. On the other 
hand, they may be used with a known mechanism to compute practical combustion 
quantities such as 

(1) Kinetic rocket or jet engine performance 

(2) Composition of engine exhaust gases 

(3) Ignition delays In combustible mixtures 

(4) Ignition lengths and nozzle performance in supersonic combustion 

Examples of some of these computations performed at NASA Lewis are given later 
In this paper. 


GENERAL THEORY 

A reacting system containing NS species Is assumed to have NR Inde- 
pendent chemical reactions proceeding simultaneously. The general set of 
reactions Is written 

NS 


1=1 

Using the law of mass action, we can write the forward and reverse rate of 
each reaction as 


NS 



where p Is the mass density of the mixture and the concentration of 
species Si, (moles per unit mass of mixture). The rate constant for the for- 
ward reaction kj Is written In the modified Arrhenius form 

n -E./RT 

kj = AjT J e 3 (3) 




j=l ,2, . . . , NR 
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The reverse rate constant k_j 
reverslbl lity: 


In these last equations, 

Ej = activation energy for the 
T = temperature, K 
R = universal gas constant 
(K e q)j = equilibrium constant for reaction j 

For any reacting system, batch or flow, the appropriate continuity equa- 
tions are differentiated to give a set of differential equations for the vari- 
ation of the concentrations and temperature with time. For a flow reaction, 
differential equations for density and velocity will also be obtained. The 
resulting system of equations Is solved numerically In combination with the 
Ideal gas equation of state: 


Is obtained from the law of microscopic 


k _j * K eq*j 


( 4 ) 


jth reaction, energy/mol 


P 

p = NS 

XT £ «, 
1=1 


( 5 ) 


where p Is the pressure of the mixture. 

For the case of constant pressure and adiabatic batch reaction, the equa- 
tions are as follows: 


d „ w 

dt - "7 U1 ' 2 " S 

where W-j , the molar rate of formation of species , Is given by 

NR 

W 1 ■ L ( *ij ~ ”13 )(R 3 " R -3’ 

3=1 

and 


dl 

dt “ ” NS 


NS 

w i h i 


E cr. C , 

i-i 1 Pr1 


( 6 ) 


( 7 ) 


( 8 ) 
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where 


= molar enthalpy of species 
c p> i = molar heat capacity of species 

The standard Runge-Kutta or other explicit numerical techniques cannot be used 
to solve these equations. Although the differential equations are Inherently 
stable, they have widely different time constants for relaxation to their 
final equilibrium solution. Therefore, the step sizes required for solution 
stability are prohibitively small. These systems of equations are called 
stiff, and the problem of solving them accurately has been studied for many 
years. Several new Implicit Integration methods have been developed for solv- 
ing these stiff equations. Although they have been quite successful In 
Improving the accuracy and efficiency of the solutions, there Is still a need 
for additional development of better methods. Problems of accuracy and of 
efficient step-size selection can still arise. Methods need to be developed 
for automatic selection of the optimum Integration parameters. Therefore, 
research Is actively proceeding at many Institutions to develop more efficient 
Integration methods (refs. 1 to 4). 

There Is another kinetic combustion model of Interest to practical engine 
designers - namely, highly backmlxed reacting flow In what Is usually called 
the well-stirred reactor. We assume the limiting condition of zero-dimensional 
flow or Instantaneous backmlxlng of the reacted gases with the cold, unreacted 
gases. Although this model Is an oversimplification of highly turbulent 
reacting flow, It Is a very useful first-order approximation for some practi- 
cal reactlng-flow systems. 

The process Is a constant pressure combustion with mass flow rate m 
through a reactor of constant volume v. The average residence time In the 
reactor Is 


For each species we can write a continuity relation which says the fol- 
lowing: The difference between the species molar flow rate Into and out of 

the reactor Is equal to Its rate of formation (or destruction) by chemical 
reaction In the reactor. These continuity equations are 

NR 

® <°1 - «i) + XI (v ij " v 1j )(R j ' R -j } 1=1 ’ 2 NS (10) 

The following energy conservation equation can also be written (If we assume 
that the process Is adiabatic): 

NS 

E (o^ - o*h*) =0 (11) 

In these equations * Indicates the unreacted gas mixture. Equations 
(10) and (11) are a set of nonlinear, algebraic equations for and the 
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reactor temperature. They are solved by the Newton-Raphson Iterative proce- 
dure. Logarithmic Increments of the variables are used to avoid numerical 
problems (ref. 5). 


NASA GENERAL CHEMICAL-KINETICS CODE 

Many computer codes have been published which perform different kinds of 
chemical-kinetic computations. At NASA Lewis we published a general chemical- 
kinetics code, GCKP84 (ref. 6). It performs a wide variety of chemical- 
kinetics computations with convenience and efficiency. 

It Is designed to perform chemical-kinetics computations for several 
reaction models. Including the following: 

(1) General reaction In either a batch system or a one-dimensional, 
frlctlonless plug flow 

(2) Combustion reaction In a well-stirred reactor (highly backmlxed flow) 

(3) Reaction behind a shockwave with boundary layer corrections 

(4) Ignition processes In either a batch or flow system 

(5) Nozzle expansion reactions 

For each of these models the following general features of the code apply: 

(1) Any chemical system may be used for which reaction rate constant data 
and species thermodynamic data are known. 

(2) The process may be adiabatic, or the heat transfer between the 
reaction and Its environment may be considered. 

(3) A new efficient Integration technique Is employed (ref. 7). 

(4) Any chemical reaction of the form aA + bB ^dD + eE may be used, 
Including photochemical and Ionic reactions. 

(5) Simplified Input for combustion reactions may be used. 

The code contains an option to compute rocket performance parameters for 
nozzle flow and can be conveniently used to compute Ignition lengths In super- 
sonic flow. An option Is also provided to perform a well-stirred reactor 
computation and, then. Immediately expand the products through an assigned 
nozzle profile. 


EXAMPLES OF CHEMICAL KINETIC COMPUTATIONS WITH GCKP84 
Hydrogen-Oxygen Batch Reaction at Constant Temperature 

Measurements of hydrogen peroxide formation In the constant-volume. Iso- 
thermal reaction of hydrogen and oxygen at 500 °C were reported many years ago 
by Baldwin (ref. 8). The Important reactions In the hydrogen-oxygen mechanism 
are as follows: 


H 2 + 0 2 = OH + OH 
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H 2 + OH = H 2 0 + H 
H + 0 2 = OH + 0 
0 + H 2 = OH + H 
*H + 0 2 + M = H0 2 + M 
H + H0 2 = OH + OH 
*H 2 0 2 + M = OH + OH + M 
*ho 2 + H0 2 = H 2 0 2 + o 2 
*H 2 + H0 2 = H 2 0 2 + H 

The rate constants for many of these reactions have been measured fairly 
accurately, but some of the reactions Involving H 2 0 2 and H0 2 still have a sig- 
nificant uncertainty. The reactions Important In determining H 2 0 2 concentra- 
tion are Indicated by an asterisk. GCKP84 was used to compute the H 2 0 2 versus 
time profile. After making several variations of the rate constants for the 
two most uncertain and Important reactions (H 2 + H0 2 H 2 0 2 + H and H0 2 + 

H0 2 ^ H 2 0 2 + 0 2 ), we obtained the agreement between the experimental and com- 
puted results shown In figure 1. This example Illustrates the use of kinetics 
computations In obtaining better rate-constant values for one or two key reac- 
tions when all other reactions In a mechanism are fairly well known. Of 
course, we cannot say that we have uniquely determined both reaction rate con- 
stants since we have two adjustable parameters. But we have put limits on 
their uncertainty. Of course, this type of exercise can only be done when the 
experimental variable being matched by computation Is sensitive to small vari- 
ations In the uncertain rate constants. This brings up the following general 
question: For a given change In any rate constant of the reaction mechanism, 

how much does each computed variable change? This Is the relatively new 
field of study In chemical kinetics called sensitivity analysis. 

In practice, even through a chemical mechanism may contain more than 100 
reactions, only the rate constants of 15 or 20 may significantly affect the 
computed results when their values are changed. Recently, methods have been 
developed to systematically compute these effects In the form of sensitivity 
coefficients (refs. 9 to 11). This computation Is performed along with the 
chemical-kinetics computation. The results of the latter computation are used 
as Input to solve a second set of differential equations whose unknowns are 
sensitivity coefficients of the form 


T 



k 0 k 
8. m 


Hm 


i c i 

ak„ 


2 

a c. 


ak ak 

Si m 


( 12 ) 

(13) 


The first-order coefficient effectively gives the percentage change In 

c^ ( 1 th concentration variable) for a given percentage change In the rate 
constant of reaction 8. The second-order coefficient Tia.m gives the effect 
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on of changes In both k a and k m . It is now recognized that a sensi- 

tivity analysis must be performed for any complex reaction In order to pinpoint 
the Important reaction paths and obtain a good understanding of the reaction 
mechanism. 


Computation of Rocket Performance 

One of the first uses of chemical kinetic computations was In the compu- 
tation of rocket performance using a one-dimensional flow model. This per- 
formance Is measured by the velocity of the hot gas at the nozzle exit point. 

As the gas expands and cools. Its velocity (i.e., kinetic energy) Is kept as 
high as possible by chemical recombination reactions which occur In the noz- 
zle. These reactions convert high potential-energy atoms and radicals Into 
low-energy stable molecules. Maximum performance Is obtained If these reac- 
tions maintain equilibrium conditions in the very short nozzle-residence time. 
But this Is often not the case since the finite rates of the reactions are not 
fast enough. By knowing the rate constants for the Important recombination 
processes one can compute kinetic performance that more accurately and realis- 
tically reflects the experimentally measured performance. An example of this 
is shown In figure 2 for the oxygen difluoride - dlborane system (ref. 12). 

The three curves of performance as a function of oxidant/fuel ratio (fig. 2(b)) 
show the maximum and minimum theoretical performance (assuming equilibrium and 
frozen conditions) as well as the klnetlcally limited performance. The latter 
agrees rather well with actual measurements corrected for various losses. The 
reaction mechanism used Is shown in table I. This system Is Interesting 
because only one reaction rate constant has to be accurately known In order to 
compute the kinetic curve. This is the recombination of hydrogen atoms to 
give molecular hydrogen. This fact was determined by a simple sensitivity 
analysis. Computed performance is unchanged for wide variations of the other 
rate constants. The reason for this situation can be seen from the heat 
release data In table I. The hydrogen atom recombination accounts for well 
over 50 percent of the heat release in the process, and this, of course, con- 
trols performance. 


Two-Stage Well-Stirred Reactor 

The practical use of the well- stirred reactor model of GCKP84 is Illus- 
trated by Its simulation of a two-stage turbulent combustor. A theoretical 
and experimental study of two-stage, rich- lean combustion as a method of con- 
trolling nitrogen oxide emissions was performed at Lewis. This technique was 
suggested as a way of reducing the conversion of the organic nitrogen compounds 
in hydrocarbons (fuel-bound nitrogen) to nitric oxide during combustion. The 
objective of the work was to study the effect of operating conditions on emis- 
sions of both nitrogen oxides (NO x ) and carbon monoxide (CO) (ref. 13). 

Both the nitrogen and the hydrogen content of the fuel were changed to simu- 
late coal-derived syncrude fuels. These fuels contain high percentages of 
aromatic hydrocarbons and therefore have a lower hydrogen-to-carbon ratio than 
today's petroleum-based fuels. Experiments were performed in the two-stage 
flame-tube apparatus shown In figure 3. The fuels were mixtures of propane, 
toluene, and pyridine blended to give a range of hydrogen and nitrogen con- 
tents. The primary-zone equivalence ratio was varied from 0.7 to about 1.8. 

The products of first-stage combustion were rapidly diluted with air Injection 
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and burned at a final equivalence ratio of 0.5. The final N0 X and CO concen- 
trations were measured. The two-stage flame tube was modeled analytically by 
a two-stage stlrred-reactor computation closely simulating the experimental 
fuels. Instantaneous mixing of the secondary dilution air was assumed. A 
detailed chemical mechanism was used for propane and toluene oxidation and for 
N0 X formation. The fuel-bound nitrogen was Input as nitrogen atoms and all 
mixing and heat-transfer effects were neglected. This very simplified model 
of a quite complex process was able to predict most of the observed trends of 
N0 X and CO formation with variation of nitrogen and hydrogen concentrations 
In the fuel. The model's prediction of other trends, which were not measured 
experimentally, gave additional Information about the advantages and disadvan- 
tages of using this rich-lean combustion for emissions control. A comparison 
of some experimental and computed results for N0 X and CO concentration as a 
function of primary equivalence ratio Is given In figures 4 and 5. The simple 
model qualitatively predicts the observed trends. It can be seen that using a 
rich primary equivalence ratio does reduce N0 X formation very significantly. 
However, at the same time, this technique Increases the formation of CO, which 
Is also undesirable. Thus, one set of conditions cannot minimize both pollut- 
ants, and tradeoffs will have to be made to obtain desired emissions levels 
for any practical situation. This work has shown that a simple kinetic model 
can sometimes be used to qualitatively explain the observed trends In a highly 
complicated combustion system. 


Ignition Delays In Benzene-Oxygen-Argon Mixtures 

Our current work Involves the study of the mechanism of hydrocarbon oxida- 
tion. We are studying the aliphatic hydrocarbons by measuring Ignition delays 
and concentration profiles during the shock-heated oxidation of hydrocarbon- 
oxygen-argon mixtures. Then a detailed chemical mechanism Is formulated and 
used In detailed kinetics computations to match the observed experimental data 
over a wide range of Initial composition, temperature, and pressure. We are 
also studying the simplest aromatic, benzene. A partial mechanism for benzene 
oxidation Is as follows: 


c 6 h 6 + °2 = c 6 h 5 + h0 2 
C 6 H 6 = C 6 H 5 + H 

c 6 h 6 = c 4 h 4 + c 2 h 2 
H + 0 2 = OH + 0 

H + 0 2 + M = H0 2 ♦ M 
0 + C 6 H 6 = C 6 H 5 + OH 
OH + C 6 H 6 = C 6 H 5 + H 2 0 
C 6 H 5 + ° 2 = C 6 H 5 0 ♦ 0 
C 6 H 5 0 + M = C 5 H 5 + CO + M 
CO + OH = C0 2 + H 
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The experimental program Is described by T.A. Brabb's paper presented at this 
symposium. I attempted to match experimental Ignition delay data for benzene- 
oxygen-argon mixtures that were taken by another Investigator at Lewis. The 
reaction zone behind the experimental reflected shock Is approximated by a 
constant volume batch reaction. The theoretically computed temperature and 
pressure (for nonreacting conditions) behind the shock are used as starting 
conditions. A typical, computed profile of pressure as a function of time Is 
shown In figure 6. The Ignition delay t Is determined from the first sig- 
nificant pressure rise; this Is similar to determining the experimental t 
value from the measured pressure trace. The mechanism we used Includes over 
100 reactions Involving benzene and Its degradation products, acetylene, phe- 
nol, phenylacetylene, ethylene, and methane. It Includes all the reactions of 
the hydrogen-oxygen system as well. Values from the literature were used for 
all reaction rate constants except the quite uncertain ones of the benzene, 
phenol, and phenylacetylene reactions. These were varied within reasonable 
limits to get the best overall agreement between computed and experimental 
t values. Starting mixture equivalence ratios ranged from 0.5 to 2, tempera- 
tures ranged from about 1200 to 1700 K, and Initial pressures ranged from 2 to 
6 atm. A simple sensitivity study made by changing rate constants one at a 
time has shown that the x values are most sensitive to the rate constant for 
the C5H5 + O2 reaction, which Is quite uncertain. Other rate constants which 
are uncertain and have a significant effect are those for C6H6 * C4H4 + C2H2 
and C&H5 C4H3 + C2H2. In figure 7 we show some comparisons between com- 
puted and experimental Ignition delays. We have used only the experimental 
data points which are considered the most accurate for the mechanism matching. 
These are the data with x values of 100 msec or greater. Shorter Ignition 
delays may have significant error due to the nonuniformity of the reaction 
mixture and the heating effects. The computed line for log x versus 1/T 
comes within ±50 percent for all but one point which differs by about a factor 
of two. While this mechanism does fairly well. It Is certainly not complete 
and has two significant problems. First, It tends to predict delays longer 
than the experimental Ignition delays for the lean mixtures and shorter than 
the experimental delays for the rich mixture. Second, It does not predict the 
experimentally observed effect of dilution with argon at a constant equiv- 
alence ratio of 1.0. Therefore, more work needs to be done on this oxidation 
mechanism. A full understanding of the oxidation will probably not be achieved 
until further experimental data are obtained on concentration profiles of a key 
species, such as the phenyl radical or carbon dioxide. Work along these lines 
Is proceeding at Lewis and other laboratories. 
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TABLE I. - 0F 2 -B 2 H 6 MECHANISM 
[Equivalence ratio = 1.63.] 


Reaction 

Net reaction 
rate, 
percent 

Heat of 
reaction, 
AH, 

kcal/mole 

H + H + Mt±H 2 + M 

32.3 

-104 

B0 + HF ^ BOF + H 

26.8 

-33 

h + f + m^hf + m 

5.1 

-136 

H + OH + M ^ H 2 0 + M 

5.2 

-119 

BF + OH ^ BOF + H 

3.9 

-58 

F + H 2 HF + H 

4.5 

-32 

H 2 + OH ^H20 + H 

9.7 

-15 

BF + 0 BO + F 

9.3 

8 

H + BF 2 r^HF + BF 

2.4 

-25 

BO + F + M ^ BOF + M 

.3 

-169 

BF + 0 + M ^ BOF + M 

.2 

-160 

O + H + M^OH + M 

.3 

-102 

Total 

100.0 



Energy release 
rate, 
percent 


100.0 


H 2 O 2 CONCENTRATION 







EXPERIMENTAL 


COMPUTED 


TIME, sec 

Figure 1. - Computed and experimental results for hydrogen -oxygen reaction. 
Temperature, T, 500 °C. 














(a) Rocket engine conditions for OF 2 -BoH 6 propellant Equivalence ratio, 
1.63. 




2.0 Z4 2. 8 3.2 3.6 4.0 

OXIDANT-FUEL RATIO 


(b) Performance of OFo-B^ rocket propellant Chamber 
pressure, 689 kPa 000 psia). 

Figure 2. - Oxygen difluoride-diborane (OF 2 -B 2 H 5 ) rocket 
combustion. 
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4000 p- 


EXPERIMENTAL 

DATA 



11% HYDROGEN 
9% HYDROGEN 
SRC COAL 
SYNCRUDE BLEND, 
8.7% HYDROGEN 
0.95% FBN 

COMPUTED RESULTS 


1.0 1.2 1.4 1.6 

PRIMARY EQUIVALENCE RATIO, <p p 

Figure 4. - Effect of hydrogen content on NO x emissions. Fuel-bound 
nitrogen, 1.0 percent; secondary equivalence ratio, 0.5; secondary 
residence time, = 2 msec. 


PERCENT 



PRIMARY EQUIVALENCE RATIO, <p p 

Figure 5. - Effect of percent hydrogen on carbon monoxide 
emissions for hydrocarbon combustion. Fuel-bound ni- 
trogen, 0 percent; secondary equivalence ratio, 0. 5 ; 
secondary residence time, s 2 msec. 
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4.0 

3.6 


3.2 



TIME, Msec 


Figure 6. - Computed pressure as function of time for benzene-oxygen-argon shock ignition. Initial 
temperature, 1366 K; equivalence ratio, <p, 2.0. 


EQUIVALENCE 

RATIO, 

<P 


□ 0.5 

2000 h A 1.0 



Figure 7. - Ignition delay data for benzene-oxygen-argon mixtures. 


